Superfluidity and collective modes in Rashba spin-orbit coupled Fermi gases 
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We present a theoretical study of the bulk superfluid properties and the collective modes in two-component 
Fermi gases with .s-wave attractive interaction and synthetic Rashba spin-orbit coupling. Both the systems in 
three and two spatial dimensions are investigated. The general effective action for arbitrary spin-orbit coupling 
and Zeeman field are obtained from the functional path integral formalism. At zero Zeeman field, the system 
undergoes a crossover from the ordinary fermionic superfluid to the Bose-Einstein condensation of a novel 
type of bound molecules (referred to as rashbons) when the spin-orbit coupling strength is tuned from weak 
to strong. We show that the bulk superfluid properties and the properties of the collective modes manifest this 
crossover. At large spin-orbit coupling, the superfluid density and the velocity of the Goldstone sound mode 
become universal, that is, independent of the strength of the .s-wave attraction. We also determine the two-body 
interactions among the rashbons. At nonzero Zeeman field, the system undergoes quantum phase transitions 
to some exotic superfluid phases which are topologically nontrivial. For the two-dimensional system, we show 
that the quantum phase transition is related to the gapless fermionic spectrum which causes infrared singularities 
and the phase transition is of third order. The superfluid density and the sound velocity behave oppositely in the 
normal and the topological superfluid phases: They are suppressed by the Zeeman field in the normal superfluid 
phase, but get enhanced in the topological superfluid phase. The sound velocity also shows non-analytical 
behavior across the phase transition. Around the phase transition, the massive Higgs mode is softened. For 
the three-dimensional system, the singularities caused by the gapless fermionic spectrum are weakened and the 
phase transition is of higher than third order. However, the superfluid density and the sound velocity have similar 
behaviors to the two-dimensional system. 

PACS numbers: 03.75.Ss, 05.30.Fk, 67.85.Lm, 74.20.Fg 



I. INTRODUCTION 



It was proposed many years ago fljj, |2j that, by tuning the 
strength of the attractive interaction in a many-fermion sys- 
tem, one can realize a smooth crossover from the Bardeen- 
Cooper-Schrieffer (BCS) superfluidity at weak attraction to 
Bose-Einstein condensation (BEC) of difermion molecules 
at strong attraction y|-|9|]. One typical example is the dilute 
Fermi gas in three dimensions with short-range attractive in- 
teraction, where the effective range ro of the interaction is 
much smaller than the inter-particle distance characterized by 
where kp is the Fermi momentum in the absence of inter- 
action. The attraction strength is characterized by a dimen- 
sionless parameter 1 /{kpa s ) where a s is the s-wave scattering 
length of the short-range interaction. The BCS-BEC crossover 
has been confirmed in the experiments of ultracold fermionic 
atoms OUnilh where the s-wave scattering length and hence 
the parameter l/(kpa s ) was tuned by means of the Feshbach 
resonance. 

On the other hand, the effect of a nonzero Zeeman field 
h has been a longstanding problem of fermionic supercon- 
ductivity/superfluidity for several decades lfl3ll . It is gener- 
ally believed that the superfluidity is completely destroyed 
when the Zeeman field becomes large enough. The well- 
known theoretical result for s-wave weak coupling supercon- 



ductors is that, at a critical Zeeman field /?cc = Ao/ V2 (called 
Chandrasekhar-Clogston limit) where Ao is the zero tempera- 
ture gap at h = 0, a first order phase transition from the BCS 
state to the normal state occurs ITil [l5ll . Further theoretical 
studies showed that the inhomogeneous Fulde-Ferrell-Larkin- 
Ovchinnikov state may survive in a narrow window between 
hcc and /ifflo ~ 0.754A d. The Zeeman field effects in 
the BCS-BEC crossover have been experimentally studied by 
using cold fermionic atoms ifnl [l8ll . The atom numbers of 
the two lowest hyperfine states of 6 Li are adjusted to create 
a population imbalance which simulates effectively the Zee- 
man field h. The experimental results show that the fermionic 
superfluidity in the BCS-BEC crossover regime is also com- 
pletely destroyed when the Zeeman field is large enough. 

The recent experimental breakthrough in generating syn- 
thetic non-Abelian gauge field and synthetic spin-orbit cou- 
pling Il9l - l26ll has opened up the way to study the spin-orbit 
coupling effects as well as the combined spin-orbit cou pling 
and Zeeman field effects on the BCS-BEC crossover II271 - 
l34l,l36ll . In solid state systems, it was shown that the topologi- 
cally nontrivial superconducting phase appears in spin-orbit 
coupled system if the Zeeman field is turned on I37l - l44tl . 
For neutral atoms, the spin-orbit coupling can be gen erated 
through a synthetic non-Abelian gauge potential A I22I1 . The 
well-known Rashba spin-orbit coupling for spin- 1/2 fermions 
can be generated via a 2D synthetic vector potential 1I24I.I25I1 



A = —Aticr ± — —Afi(cr x e x + cr v e v ), 



(1) 
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where a ± = a x e x + a v e v for any vector a. The single-particle 
Hamiltonian for a fermion moving in the synthetic gauge field 
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is given by "Hq — (p - A) 2 /(2m) where and p = —ihV is the 
momentum operator. In this paper we use the natural units 
h — = tn — 1 for convenience. 

For the 2D synthetic vector potential A given in (Q]), the 
single-particle Hamiltonian can be reduced to 



•Ho = — + Ao- ± ■ p ± , 



(2) 



where an irrelevant constant A 2 /2 has been omitted. The 
spin-dependent term Axr± • p_i_ can be mapped to the standard 
Rashiba spin-orbit coupling A(cr x p y - cr y p x ) by a spin rotation 
o~x — » °"v an d °"v - > _cr .v The gauge field strength i char- 
acterizes the strength of the spin-orbit coupling, which can be 
tuned from weak to strong in cold atom experiments. Since 
the final physical results depend only on the A 2 , we set A > 
in this paper without loss of generality. For many-fermion 
system, the spin-orbit coupling strength can be characterized 
by the dimensionless ratio A/kp. While for solid state systems 
this ratio is v ery small, it can reach the order 0( 1 ) in cold atom 
systems 12(3. 1210 . Therefore, cold fermionic atoms provide the 
way to study the fermionic superfluidity from weak to strong 
spin-orbit coupling. 

Motivated by the experimental progress of realizing spin- 
orbit coupled atomic Fermi gases, the fermionic superfluid- 
ity with spin-orbit coupling has been extensively studied l27l - 
[34l [36h . It was shown that, in the presence of Rashba spin- 
orbit coupling, the two-body bound state exists even for a s < 
where the bound state does not exist when A — 0. With in- 
creased A, the binding energy is generally enhanced J45ll . The 
bound state also possesses a non-trivial effective mass which 
is generally larger than twice of the fermion mass l2& [29T,l32Tl . 
Such a novel bound state is referred to as rashbons in the lit- 
eratures l46ll . For many -body system, it has been proposed 
that a spin-orbit-coupled Fermi gas can undergo a smooth 
crossover from the ordinary fermionic superfluid to the Bose- 
Einstein condensation of rashbons even for weak attraction if 
A/kp is tuned from small to large values |27U33ll . On the other 
hand, if a nonzero Zeeman field h is turned on, some topolog- 
ically nontrivial superfluid phases can emerge Sii. 

In this paper, we study the bulk superfluid properties and 
the collective modes in Rashba spin-orbit-coupled Fermi su- 
perfluids. We mainly consider two aspects: (1) the bulk super- 
fluid properties and the collective modes from weak to strong 
spin-orbit coupling at zero Zeeman field, which manifest the 
crossover from ordinary Fermi superfluid to the Bose-Einstein 
condensation of rashbons and (2) the quantum phase transi- 
tions from normal superfluid phase to topologically non-trivial 
superfluid phases in the presence of nonzero Zeeman field and 
their effects on the bulk superfluid properties and the collec- 
tive modes. 

This paper is organized as follows. In Sec. II, we derive the 
general effective action for the superfluid ground state and the 
collective modes with arbitrary spin-orbit coupling and Zee- 
man field by using the functional path integral method. In 
Sec. Ill and IV, we study the systems with zero Zeeman field 
in three and two spatial dimensions, respectively. The sys- 
tems with nonzero Zeeman fields are studied in Sec. V. We 
summarize in Section VI. 



II. GENERAL FORMALISM 

We consider a homogeneous Fermi gas with a short-range 
i-wave attractive interaction in the spin-singlet channel. For 
cold atom experiments, the attractive strength can be tuned 
from weak to strong [80 . In the dilute limit where the effec- 
tive range ro is much smaller than the characteristic length 
scales of the system, that is, ro <k k F l , a s , A~ l , the interac- 



tion Hamiltonian can be modeled by a contact interaction (4 
The many-body Hamiltonian of the system can be written as 



H — Ho + Hz + H n 



(3) 



where 



H 
H z 



, (r)[^-+Ao- ± -p ± -ij\^(rl 



= j d 3 rif/\ 

= -h J cPrt/j\r)cr z t/j(r), 
H int = -U J d 3 r ^(r)^(r)^(r)^ T (r). (4) 

Here tf/(r) = [^j(r), i/^(r)] T represents the two-component 
fermion fields, /u is the chemical potential, and h is the Zee- 
man magnetic field. We set h > in this paper without loss of 
generality. The contact coupling U > denotes the attractive 
s-wave interaction between unlike spins. 

In the functional path integral formalism, the partition func- 
tion of the system is 



Z= Dt//Dt//exp{-S[i(r,i(r]}, 



where 



rP r rP 
S[i/r, tfr] = I dr I d 3 rt//d T i// + I drH(t//, 
Jo J Jo 



(5) 



<A). (6) 



Here /3 — l/T and H(tf/, iff) is obtained by replacing the field 
operators iff* and if/ with the Grassmann variables \/i and if/, re- 
spectively. To decouple the interaction term we introduce the 
auxiliary complex pairing field <D(x) = -Uif/i(x)if/f(x) [x = 
(t, r)] and apply the Hubbard-Stratonovich transformation. 
Using the Nambu-Gor'kov representation 



(*$(,))• *C*H*« -^Mto-,). (7) 



we express the partition function as 

Z= j OTSTOO^XD'expJ-^t^^^d) 1 ]}, 



(8) 



where 



JPF,*,<M>1] = I j dx\®(x)\ 2 
~\f dx f dx"y(x)G- y (x,xy¥(x). 



(9) 
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The inverse single-particle Green's function G l (x,x > ) is 
given by 



G _1 (x,x') = 



G- l (x) <D(jc) 

o'Xx) G: 1 ^) 



5(x - .x'), 



(10) 



where 

G^(jc) = -d T + hcr z + (p 2 /2 + A<r ± ■ p±-/j). (11) 
Integrating out the fermion fields, we obtain 

Z= j D<$>D& exp { -S eff [<l>, O f ]}, (12) 

where the effective action reads 



dx\®(x)\ 2 - -TrlntG" 1 ^, x')]. (13) 



The effective action .SefffO, O*] cannot be evaluated pre- 
cisely. In this work, we consider mainly the zero tempera- 
ture case. Therefore, we follow the conventional approach 
to the BCS-BEC crossover problem, that is, we first consider 
the superfluid ground state which corresponds to the saddle 
point of the effective action, and then study the Gaussian fluc- 
tuations around the saddle point. The Gaussian fluctuations 
corresponds to the collective modes, including the gapless 
Goldstone mode and the massive Higgs mode. In ordinary 
fermionic superfluids, only the Goldstone mode remains at 
low energy whereas the Higgs mode is pushed up to the two- 
particle continuum. Therefore, the Higgs mode usually ap- 
pears as a broad resonance at the large characteristic energy 
scale of the system. 

In the superfluid ground state, the pairing field <D(x) ac- 
quires a nonzero expectation value (<t>(x)) = A, which serves 
as the order parameter of the superfluidity. Due to the U( 1 ) 
symmetry, we can set A to be real without loss of generality. 
Then we separate the pairing field as <D(x) = A + <p(x), where 
<p(x) is the fluctuation around the mean field. The effective ac- 
tion iSefftO, O' ] can be expanded in powers of the fluctuation 
<p(x), that is, 



SefftO, O 1 ] = Sf^(A) + S^[cf>, f] + 



(21, 



(14) 



where 5^ (A) = tS e ff[A, A] is the saddle-point or mean-field 
effective action with the superfluid order parameter A deter- 
mined by the saddle point condition OS 2 /dA = 0. Note that 
under the saddle point condition the linear terms in <p and ' 
in Eq. ( TT4T > vanish automatically. 



A. Saddle point: mean-field approximation 

The mean-field effective action or the thermodynamic po- 
tential Q. can be expressed as 



£1 



f3V 



= —--J]\ndet[g- l (K)], (15) 



where the inverse fermion Green's function reads 

6 (K) -\ A @Z\K) 
and Q~^(K) is given by 

@l\K) = ia>„ + hcr : + (£ k + A(r ± ■ kj. 



(16) 



(17) 



The dispersion £k is defined as ^ = ek - ft with ek = k 2 /2. In 
this paper K = (io) n , k) denotes the energy and the momentum 
of fermions with lo„ = {2n+\)nT {n integer) being the fermion 
Matsubara frequency. We use the notation YiK - | Zn Zk 

with 2k = fd 3 k/(2n) 3 for 3D system. 

The determinant of the inverse fermion propagator, 
&et[Q~ x {K)], can be evaluated as 

] [(ia>„ + shf -El- A 2 k 2 x ] - 4^ 2 ki(^ 2 - h 2 ) 

s=± 

= [(icj n ) 2 - (£+) 2 ] [(iu„) 2 - (E k ) 2 ] , (18) 



where Ek = -J^ + A 2 and the quasiparticle excitation spectra 
are given by 



2&. (19) 
Here the quantities rjk and & are defined as rjk = -^/l 2 k 2 + h 2 



and = yj^k^k + ^ 2 ^ 2 - Completing the Matsubara fre- 
quency sum we obtain the explicit form of the mean-field ef- 
fective action 



A 2 

Q = — 

U 



zz 

k a=± 



pa 
£ k 



-^ + rin(i + - 



(20) 



where = £ k ± Ik- Here the term \ 2k ^ = 2k ft is 
added to recover the correct limit for A — 0. The integral over 
the fermion momentum k is divergent and the contact cou- 
pling U needs to be regularized. For a short-range interaction 
potential with its s-wave scattering length a s , it is natural to 
regularize U by the two-body T-matrix in the absence of SOC. 
We have 



1 

U 



1 

47Tfl t 



(21) 



The superfluid order parameter A should satisfy the saddle 
point condition <90/<9A = 0, or the so-called gap equation 



1 

U 



ZZ' 1+ « L 

k a=± 



- 2f(E«) 



(22) 



where f{E) = l/(e E ^ T + 1) is the Fermi-Dirac distribution. 
Meanwhile, if the total fermion density n is imposed, the 
chemical potential [i should be determined by the number 
equation -d£l/dfi = n, that is, 



zz 

k a=± 



1 f 1+ 

1 + a — 

2 { ft 



1 - 2/(£?) 



2E» 



(23) 
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In general, A and /j. are obtained by solving the gap and 
number equations simultaneously. As a convention, we de- 
fine the Fermi momentum kp through the noninteracting form 
n = kp/{3n 2 ), and the Fermi energy is given by ep = k 2 /2. 

In the Nambu-Gorkov space, the fermion propagator 0(K) 
takes the form 



where L ± (K) are given by 



L ± (K) = ia)„ - ha z ± (& - Aa ± • kj. 



(26) 



0(E) = 



0u(E) 012(E) 



021(E) 022(E) 

The matrix elements can be evaluated as 

L + (K)L-(K) - A 2 



(24) 



To evaluate the collective mode propagator, we also express 
the fermion propagator in an alternative form by using the fol- 
lowing projectors 



K(h) = = 1 ± 



1 /, A(r ± ■ k ± + hcr z 



0u(E) =0Z\E) 
022(E) =0l\E) 
0n(E) = -A 
02i(E) = -A 



[(iCO n ) 2 -(Ei) 2 W(On) 2 -(Ei) 2 y 

L^(K)L + (K) - A 2 
[(iaj n ) 2 -(E£) 2 ][(icj„) 2 -(E k ) 2 ]' 
L-(K)L + (K) - A 2 
[(ico n y--(E^) 2 ][(i^) 2 -(E k ) 2 ]' 

L+(K)L-{K) - A 2 
[(ito n ) 2 -(E£) 2 ][(iuj„) 2 -(E k ) 2 Y 



Ik 



which possess the following properties 



(27) 



(25) 



T + k (h) + P~(h) = 1 , PRhypkh) = Safi%{h). (28) 



With the help of these projectors, the fermion propagator can 
be expressed as 



0n(E) 

022(E) 

012(E) 



021(E) 



v [gghg - (^ k g ) 2 ]n t (- /i ) - A2y ^w 

£T^ + * J [(^„) 2 -(£ k + ) 2 ][(/ W „) 2 -(£ k ) 2 ] ' 

y [('^) 2 ~ (y) 2 ]PjIW ~ A 2 ^(-/Q 

[(/ w „) 2 -(£ k + ) 2 ][(/ W „) 2 -(£ k ) 2 ] ' 

[(igj„) 2 - (t k a ) 2 - A 2 ]^(/i) - 2MMon + g*)^^) 
[(i<o n ) 2 - (E+) 2 ][(i<o„) 2 - (E£ 2 ] 
y [(/^„) 2 - (g k g ) 2 - A 2 ]!P^-/i) - 2/ i (^„ - £?yPl(-h)<r z 
£i [(ico n ) 2 -(E k ) 2 ][(iw n ) 2 -(E k ) 2 ] 
y [(/^„) 2 - (^) 2 - A 2 ]^(/7) - 2h(iaj„ + ^°)^K 

2j [(^„) 2 -(^) 2 ]K^„) 2 -(^) 2 ] 

y [(/^„) 2 - (g k ") 2 - A 2 ]!P^(-/t) - 2/<^„ - g^)p-^(-ft) 
[(ioj„) 2 - (El) 2 ][(ioj n ) 2 - (E k ) 2 ] 



a=± 



(29) 



We note that the anomalous Green's function 0n(K) is not 
diagonal in the spin space for A + 0. Therefore, the spin- 
orbit coupling generates spin-triplet pairing even though the 
order parameter A have the s-wave symmetry. According to 
the Green's function relation, the spin-singlet and spin-triplet 
pairing amplitudes can be read from the diagonal and the off- 
diagonal components of 0n(E). We obtain 



<<0r T (k)^(k)> = -<^(k)^ T (k)> 



(30) 



for the spin-singlet pairing amplitudes and 
(^(k)^(k)) = -A(k x -ik y )A^— Ya- 



AEl 



<^(k)^(k)> = A(k x + ik y )A^— -Va 



ft -Ay l-2/(£«) 
& £ 4£ k 



(31) 



for the spin-triplet pairing amplitudes. In the mean field the- 
ory, the condensation density «o is half of the summation of 
all pairing amplitudes squared I48ll49ll . that is, 

"° = *Z Z KMk¥<r-(k)>l 2 - (32) 

k tr,<r'=T,J. 

It is also useful to reexpress the mean-field theory in the 
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helicity representation ||4J]. The helicity basis (i/' + ,i/'_) t is 
related to basis (if/p ipi) T by a k-dependent SU(2) transforma- 
tion. In the helicity basis, Hq + Hz is diagonal, that is 

H +H Z = Y J W + (k) + ^t(k)^_(k)] . (33) 



Therefore, the spin-orbit coupled Fermi gas can be viewed as 
a two-band system. The Zeeman field provides a band gap 
2/i at k = 0. In the presence of attraction, the mean-field 
approximation for //; nt reads 



where £(0 is defined as 

In this paper Q - (iv n ,q) denotes the energy and momentum 
of bosons with v„ = 2nnT being the boson Matsubara fre- 
quency. 

After taking the trace in Nambu-Gor'kov space, we find that 
S ~g can be written in a bilinear form 



H * - \ Z Z [^wUwlM + h,.} . 04) £^!i = 1 2 ( 0t(0 )M(0 1 *a j (40) 



o- i 8=± k 

The new momentum-dependent pair potentials A a p(k) read 



A + _(k) = -A_ + (k) = -A s (k), 
A ++ (k) = AI_(k) = -A t (k), 



(35) 



where the interband and the intraband pair potentials are given 
by 



h A(k x - ik v ) 

As(k) = —A, A,(k) = — -A. 

?7k J]k 



(36) 



Using these new pair potentials, the quasiparticle spectra 
can be expressed as 

E£ = ^(^ + |A s (k)| 2 "± + | At (k)|2. (37) 

The above expressions in the helicity basis will help us under- 
stand some results in Sec. V. 



where the inverse boson propagator M(0 is a 2 x 2 matrix, 



W M 2 i(Q) M 22 (0 



(41) 



The matrix elements of M(0 can be expressed in terms of the 
fermion propagator Q{K). We have 



M„(0 = i + i Yj Tr + 2)0 



A' 



! 22(*Q] , 
(*)], 



M 22 (0 = i + i J] Tr [0 22 (£ + 0£„ 
M 12 (0 = 1 2 Tr + Q)6n{K)} , 

L K 

M 2i (0 = X - J] Tr [^(X + 00 21 (£)] . (42) 



B. Gaussian fluctuation: collective excitations 

Then we consider the fluctuations around the mean field. 
The linear terms which are of order 0(<p) vanishes precisely 
once the saddle point condition A = Ao is imposed. The 
quadratic terms, corresponding to Gaussian fluctuations, can 
be evaluated as 

ggjftfl _ y f 10(g)! 2 

+ Z Z Tr [&{K + }. (38) 

K > 

I 



The explicit form of these functions are evaluated in Ap- 
pendix A. It is straightforward to show that Mn(iv„, q) = 
M 22 (-!'y„,q). However, we have Mi 2 (iV„,q) + M 2 i(iv„, q) 
if the spin-orbit coupling A and the Zeeman field h are both 
nonzero. Taking the analytical continuation iv„ — > a> + i0 + , 
the dispersion cj(q) of the collective mode is determined by 
the equation 

detMMq),q] = 0. (43) 

We can decompose Mn(w,q) as Mn(w,q) = M^(w,q) + 
M u (u>, q), where Mjj(o>, q) and M n (w, q) are even and odd 
functions of co, respectively. Their explicit form read 



M+( W ,q) = i + 5ZX>f(M) 



(0-E° -El 

k+p k-p 



lj + E: 



k+p k-p 



l-/(^k + p)-/(< p )l 



4 ^ x<(k, q) 



arjS=± k 



k+p k-p k+p k-p 



[f(K +P )-f(ElA- 



(44) 
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where p = q/2 for convenience, and 

afi=± k 

4 £ 2>f(M) 



1 



1 



a£=± k 



k+p k-p 
1 



co + E" + £f 

k+p k-p 

l 



o + E?. 



co -El +El 

k+p k-p 



fi-/(^ +P )-/(< P )] 
[/(^ +P )-/(< P )]- 



J k+p k-p 

On the other hand M 12(01, q) and M2i(o>, q) are even functions of co and can be expressed as 

M, 2 (w, q) = M\ 2 (co, q) + iM[ 2 {io, q), M 2 i(w, q) = M + 12 (co, q) - iM n (co, q). 
Here M j 2 (w, q) and M~ 2 (co, q) read 



(45) 



(46) 



M? 2 (a»,q) = -5 Z Z^M) 

aj8=± k 

4 ^ 2 <(k )q) 



<i./j=± k 



and 



Mr 2 (",<D = 4 Z Z^(M) 

ff,/3=± k 

4 X Km 



it,/3=+ k 



1 




1 




Hp 


co + E" +Ei 

k+p k-p 


1 




1 


co + K* - 

L k+p 


k p 


co -El +Ei 

k+p k-p - 


1 




1 


[^- £ k + p 


-E? 

Hp 


+ pa , pP 
^ + H+p + Hp 


1 




1 




Hp 


co-E" +E P . 

k+p k-p J 



[l-/(^+p)-/(^-p)] 
[m +v )-f(Ei p )] 

[i-/(^k+p)-/(<-p)] 
[/(£k+p)-/(<- P )] 



(47) 



(48) 



The explicit expressions of the functions < W" l3 (k,q) and 
1/°^(k, q) (z = 1,2,3,4) are presented in Appendix A. We 
note that <Vff (k, q) and l/^k, q) are odd functions of h, that 
is, they are proportional to hA 2 . However, the determinant of 
the matrix M is an even function of h. 

To make the results more physically transparent, we de- 
compose the complex fluctuation field <p(x) into its amplitude 
mode p(x) and phase mode 8(x), <p{x) - p(x) + iAo8(x). Con- 
verting to the variables p(x) and 9(x), we obtain 



is also called the Anderson-Bogoliubov mode. Another col- 
lective mode or the so-called Higgs mode is massive. It is 
likely heavily damped since its mass gap is generally larger 
than the two-particle continuum at q = 0. 

We are interested in the low energy behaviors of these col- 
lective modes. For this purpose, we make a small q and co 
expansion of N(0. The spin-orbit-coupling term breaks the 
0(3) rotational symmetry to a 0(2) circular rotational symme- 
try. Therefore, the expansion takes the form 



pv 



\ J] ( P*(<2) 0*02) ) N(0 ( f*® ) , (49) 



2(M| 1+ Mt 2 ) 
2A 2 (M + U - MJ 2 ) 
2A ( )M n 



A + C ± q 2 ± + Cnq? - Deo 2 



J±<1± + J\\% 
-Bco H . 



Rco 2 + 



(52) 



where the matrix N(Q) reads 
N(0 = 2 



*'Ao(M n 

-/Ao(M- + iM n ) Ag(M* 



M|,+M| 2 



■ My 



(50) 



From the expressions of M n and M^ 2 , we have M n (0, q) = 
and M n (co, 0) = 0. Therefore the amplitude and phase modes 
decouple completely at co = and q = 0. Furthermore, at the 
saddle point A = Aq we have precisely 



Mt,(0,0) = M+(0,0). 



(51) 



Here ay = a z e z for any vector a. Note that the term M 12 (o>, q) 
has no contribution to this expansion up to the order O(co 2 , q 2 ). 
We should emphasize that such an expansion is only possible 
at zero temperature, since the terms proportional to f(E£ ) - 
f(E%_ ) have the Landau theory singularity for q and co going 
to zero @I , In the following of this paper, we restrict our 
studies in the zero temperature case. 

The parameter A (J ± and /y) characterizes the stability of 
the saddle point A = Aq against the amplitude (phase) fluctu- 
ation. First, it is easy to show that 



Therefore the phase mode at q = is gapless, that is, the 
Goldstone mode. For neutral fermionic superfluids, this mode 



A=A 



(53) 
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Therefore, the stability against the amplitude fluctuation re- 
quires that A > 0. Second, the superfluid phase stiffness J ± 
(7||) is precisely proportional to the superfluid density nf 
We have 



■'J- ~~ A ' 

4m 



j _ «£_ 

" 4m 



(54) 



On the other hand, the superfluid density can be obtained from 
another equivalent definition l50l l5~lll . When the superfluid 
moves with a uniform velocity v s = v±ej_+V\\e z , the superfluid 
order parameter transforms like <1> — > <&e 2mV *' T and O* — > 
Q* e -2unv, r = 1 j n our units). The superfluid density n s is 
defined as the response of the thermodynamic potential Q to 
an infinitesimal velocity v s , i.e., 

Q(u s ) = Q(0) + + l -^\v\ + 0(wj). (55) 

Therefore, the stability against the phase fluctuation requires 
that 7 X > and J\\ > 0. Once the superfluid phase stiffness 
becomes negative, the saddle point A = Ao is unstable and 
some Fulde-Ferrell-Larkin-Ovchinnikov-like state with inho- 
mogeneous phase and/or amplitude modulation will be ener- 
getically favored ll52l - o7ll . 



As long as the stability conditions are satisfied, the disper- 
sion of the Goldstone mode at small momentum is given by 



«Xq) = ^(cf) 2 ql + (4) 2 q 2 , < 56 > 
where the transverse and longitudinal velocities are given by 



Jx 



B 2 /A+R' 



B 2 /A + R ' 



(57) 



The Higgs mode is massive and its mass gap Mh reads 



Mu = 



B 2 +AR 



DR 



(58) 



We note that the expansion ( 1521 is valid only for small fre- 
quency (l>, therefore this formula only gives the qualitative be- 
havior of the Higgs mode. In general, the Higgs mode is a res- 
onance since its spectral density arises above the two-particle 
continuum E c (q) = mink{-E k+ /2 + E, 



"k-q/2< 



Finally, we summarize the explicit expressions for the ex- 
pansion parameters in Eq. d52b . For details of the calculations, 
see Appendix A and B. First, A can be evaluated as 



5ZZ 

a=± k 



(ED 3 



h 2 Y h 4 A 2 
1 + a — + a 



+a— I 



(59) 



The expansion parameters in the frequency expansion read 

A 2 k 2 



B 



- — V V & 
4 Zj Zj (£«)3 



a± k 



l+a- 



l£E^ 

a > 



+ aY 



ft 



^(E+ k+ E k ) 2 



1 1 \ h 2 E\ 
—: + — 



1 \h 2 



1 

A 2 

R= 



zz 

a=± k 



fl + ril+la& AWjl | A 2 ,l 2 k 2 /z 2 



(£H 5 



yy 1 ^ 2 kj^ 2 y 1 

^V^k) 3 £ V( £ k+£ k ) 3 £ 1" KE k E^E-Ei) 



1 



y l . r n 

V( £ k ++£ k ) 3 a 



1 + 



K E k ) ft 



El 



k 



E k E k ) 



(60) 



The transverse phase stiffness J ± and the longitudinal phase stiffness J\\ take the form 



J± = i_L_y y Jl 

4m 1 Zj ?f? 



711 " 4m" 



2£< 

k «=± k 



k a=± 



h z Ej A 2 k 2 ± h 2 A 2 ^ 



1 + 



1 - 



2(1 



(61) 



The delta functions 6(E k ) in the expressions of J ± ,J\\ and 
A come from the zero-temperature limit of the function 
(l/47 , )sech 2 (£'^/2r) = -dfiE^/dE". The integrations over 
these delta functions vanish precisely when the excitation 
spectrum E k is fully gapped. However, at large enough 
Zeeman field h, the superfluid may become a gapless state 



where the lower excitation spectrum E k has zeros. In this 
case, these terms may have finite contributions, depending 
on whether these zeros correspond to gapless Fermi sur- 
faces ll57ll . For vanishing spin-orbit coupling, the expansion 
parameters A,B,D,R reduces to the known expressions in the 
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previous literatures 159116011 . They are given by 



where 



2 p2 
k c k 



0(£ k - h) 



E k 



6(E k - h) 



k C k 



4 V # 

k k 



(62) 



The phase stiffness becomes isotropic for A — 0. It reads 



k 2 



(63) 



1 +QjS 



/i 2 (ki - pi) - h 2 



^k+p^k-p 



(69) 



To study the two-body problem in the absence of medium 
effect, we discard the Fermi-Dirac distribution functions. The 
energy-momentum dispersion a> q of the pair excitation is de- 
fined as the solution a> + 2/j. = oj q of the equation Rer~'(w + 
ie, q) = . For h = 0, after some manipulations, the two-body 
equation becomes 



1 

U 



"Jkq 



k 6: 



kq 



i , i 2 qj sin 2 ip 



(70) 



where fikq = k 2 + q 2 /4 - w q and ^ is the angle between k ± 
and q ± . 



Therefore, for A — the zeros of E k correspond to some gap- 
less Fermi surfaces, which leads to lar ge n egative contribu- 
tions to the parameters J x , /y and A 152146011 . 



C. Two-body Problem 

To understand the behavior for the collective modes in the 
BCS-BEC crossover, it is useful to compare it with the result 
of the two-body problem in the absence of medium effect. In 
the functional path integral formalism, the two-body vertex 
function T~ l (Q) can be obtained from its coordinate represen- 
tation defined as 



r~W) = 



1 (frSeff 



PV 8<b*{x)8<$>{x') 



(64) 



4=0 



For A = 0, we have Q\2 = Qi\ - and the single-particle 
Green's function Q{K) reduces to the form 



0-{K)y 
where the diagonal elements can be expressed as 

£+(/«„, k) = V —L—r«(-h), 



(65) 



@_{iu n ,K) = Y T ^— 
Therefore, T~ l (Q) can be expressed as 



r-'(G) = i + \ Y Tr [Q + (K + Q)Q-{K)] . 



2 V 



(66) 



(67) 



The explicit expression can be evaluated as 



afi=± k ?k+p T ^k-p 



7%. «*> 



D. Two-Dimensional Case 

The above formalism applies also for the case of two spatial 
dimensions. To apply the general formalism, we only need to 
freeze the longitudinal (z) degree of freedom, that is, perform 
the following replacement: 



k 2 -> k 2 = ki, 



d 2 k 



(71) 



In the 2D case, in contrast to the anisotropic 3D case, the su- 
perfluid ground state is isotropic. A quasi-2D cold atomic gas 
can be realized by arranging a one-dimensional optical lattice 
along the axial (z) direction and a weak harmonic trapping 
potential in the radial (x - y) plane 16111 . such that atoms are 
strongly confined along the axial direction and form a series 
of pancake-shaped quasi-2D clouds. The strong anisotropy of 
the trapping potentials, namely lo z » a> ± where w z (co_i_) is 
the axial (radial) frequency, allows us to use an effective 2D 
Hamiltonian to deal with the radial degrees of freedom. 

For the 2D case, two-body bound state exists for arbitrarily 
weak attraction. The coupling constant U should be regular- 
ized in the following way l@] 



TJ £-i 



1 



U Z-J 2e k + e B ' 



(72) 



where ee is the binding energy of the two-body bound state in 
the absence of spin-orbit coupling. For convenience, we can 
define a 2D scattering length «2d by ll62tl 



4(T 2 ? 



(73) 



where y 0.577216 is the Euler's constant. For quasi- 
2D cold atoms confined by an axial trapping frequency u> z , 
the binding energy is related to the the 3D s-wave scatter- 
ing length a s by 6% — (CHa> z /7:)exp[y/2nl z /a s ] ll63ll . where 
l z = ^ZandC = 0.915. 
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III. RESULTS FOR ZERO ZEEMAN FIELD: 3D CASE 

In this section, we study the 3D system with zero Zeeman 
field (h = 0). We focus on some bulk superfluid properties 
and the collective modes from weak to strong spin-orbit cou- 
pling A. Some results for the bound state and the BCS-BEC 
crossover have been addressed in the previous literatures 1I27I - 
I29II45I1 . Here we show these results for the sake of complete- 
ness. 



Mg is always larger than 2m, and it approaches 4m for k — > 
-00. At unitary (a s — > +00) or for large spin-orbit coupling 
(A — > 00), we have k — » 0. In this case, the binding energy and 
the transverse effective mass read 



M 1 - 

77(0) = 0.439, -^- = 1.20. 

2m 



(79) 



This novel bound state is referred to as rashbons in the litera- 
tures Hi. 



A. Bound State and BCS-BEC Crossover 

First, we show that there exists two-body bound state in the 
presence of spin-orbit coupling even for negative values of a s . 
For small q, the dispersion ai q in Eq. ( l70l i can be written as 
w q = -Eb + q^/(2Mg) + qjj/(2Mg). From the imaginary part 

of the function Y~ l {to + ie, q), we conclude that the bound state 
exists if the binding energy E% > A 2 . The binding energy Eq 
is determined by the equation 



Ana. ^ 



k 2 +£ H 



(k 2 +£ B ) 2 -4/I 2 ki 



1 

k 2 



(74) 



Using the condition Eb > A 2 we find that the solution can be 
expressed as 



Eb , , , 

= 1 + m, 

where k = l/(Aa s ). Completing the integral we obtain 



r — — 1 yrr^o + i 

VI + 77W - - In = k. 

2 yi+77(K)' 



1 



(75) 



(76) 



Solving this equation, we find that a positive solution for 77 al- 
ways exists for arbitrary k. We have the exponential behavior 
tj(k) — > e 2K for k — > -00 and t]{k) — > k 2 - 1 for k — > +00. 
Therefore, bound state can form in the presence of spin-orbit 
coupling even for negative values of a s . Meanwhile, expand- 
ing Eq. ( TTOl i to the order 0(q 2 ), we obtain M fi = 2m for 
arbitrary k and the equation for the transverse effective mass 



M x 



Mir I 



-z- 



(k 2 + £ B ) 2 + 4i 2 ki 
/ir [(k 2 + £ B ) 2 -4^ 2 ki] 2 

8^ 4 ki 



E (k 2 + £ B )[(k 2 + £B) 2 -4i 2 k 2 J 2 
Completing the integral, we get 

nix) , *1(k) 1 



2m 1 

K = "2 



In 



1 + tj(k) 1 + tj(k) 1 + t)(k) 



(77) 



(78) 



The K-dependence of the binding energy and the transverse 
effective mass can be numerically obtained. The results are 
shown in Fig. Q] We find that the transverse effective mass 




FIG. 1 : The binding energy £b (we show the dimensionless quantity 
t}) and the transverse effective mass Mg (divided by 2m) as functions 
of l/(Aa s ). 



Then we turn to the superfluid medium. For zero Zeeman 
field h = 0, the single particle excitation spectra reduce to 
E^ = -y/(£k + Ak ± ) 2 + A 2 . The fermion Green's function takes 



a simple form 



-022(-z'o>„,k) = 2^ 



i0J n + ^ 



(ioj„) 2 



(E?) 2 



P£(0), 



g n (ia>„,\i) = £?2i0'«„,k) = ^ 



-A 



(iu>„) 2 



(£ ^^k(0).(80) 



where ^ = £t + and the projectors become ^(0) = 
I (1 ± cr x ■ k ± /k ± ). Since the total density n = krl/(3n 2 ) is 
fixed, the system can be characterized by two dimensionless 
parameters l/(kpa s ) and A/kp. The order parameter A and the 
chemical potential ji can be determined in units of the Fermi 
energy ep = k 2 /2 from the gap and number equations. For 
h = 0, they becomes 

1 

47ra., 

and 

1 ^ ^ 

2 



Z 

k 



y l 

^4V(^k + a^±) 2 + A 2 



1 

k 2 



£ k + aAk ± 



(81) 



(82) 



Applying the transformation k ± 
can be written in another form 



V(& + aAk ± ) 2 + A 2 

k ± +A, the above equations 



1 



4na s 



A 
An 2 



dk 7 



dk ± 



^ k -A 2 /2) 2 + A 2 



yl-L-L 

V 2£ k k 2 



(83) 
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and 



A /"*°° r*d 
n = —r I dk- I dk ± 
4n- Jo Jo 



&-^ 2 /2 



Z('-f)- 



(84) 



The integrations over k : can be analytically carried out with 
the help of the elliptic functions, which helps us to obtain so- 
lutions with high accuracy. 

The numerical results for A/e F and ^/e F are shown in Fig. 
[2] . For large spin-orbit coupling, the chemical potential be- 
comes negative and A «C \/i\. Therefore, the gap and number 
equations can be expanded in powers of A 2 /|/i| 2 . To the lead- 
ing order, the gap equation can be approximated as 



- — = 7 



\t 2 -2/j 



1 



(k 2 - 2//) 2 - 4A 2 k 2 ± k 2 



which gives 



Eb 
' 2 



Then the number equation becomes 

(k 2 + E B ) 2 + 4A 2 k 2 ± 



« = 2A 2 ^ 



[(k 2 + £ B ) 2 -4^ 2 ki] 2 



(85) 



(86) 



(87) 



which yields 



A 



16t,(k) 



3tt ^1 + t](k) k F 



For large A/kp, we have k ^ 0. Using the result 77(0) 
we obtain the following asymptotical behavior 



0.788 A E £*-1.44 

€p V €F 



(88) 



0.439, 



(89) 



The above result indicates that for A/kp — > 00, the properties 
of the system become independent of the interaction parame- 
ter l/(k F a s ). 

Beyond the leading order of A 2 /\fj.\ 2 , the chemical poten- 
tial at large A/kp can be expressed as fi = -Eg/2 + fj. B /2, 
where fi B <s E b can be referred to as the chemical potential of 
the rashbons. To determine this chemical potential as well as 
the interactions among the rashbons, we construct the Gross- 
Pitaevskii free energy of the rashbon condensate. To this end, 
we derive the Ginzburg-Landau free energy functional 



r GL [A(x)] 



= J dx A t (jc) ^1 



|fl|l-^Vi-^|V 2 ) 



+ C |A(*)i 2 + -diAoor 



(90) 



according to the fact A <& \/u\. The coefficients a, b ± , b\\ can be 
obtained from the two-body vertex function T~ l (Q) and c,d 



1 - 

< 



0.4 

0.2 






1.5 




FIG. 2: The pairing gap A (divided by ep) and the quantity /j + £b/2 
(divided by e F ) as functions of A/k F for various values of \l{k 7 a s ). 
The dashed lines corresponds to the analytical results J88b and l !94l l 
for large spin-orbit coupling. 



can be obtained from the ground state energy Q(A). To the 
leading order of fi B /E B , these coefficients can be evaluated as 



Z 



(k 2 + E B ) 2 + 4A 2 kl 1 VI + r^K) 



[(k 2 + E B ) 2 - 4A 2 k 2 J 2 A Stu](k) 
a, b\\ = -a, c = -a/j B , 



2K 



2M 



B 



,=22 

k 

1 



(k 2 + £ B )[(k 2 + E B ) 2 + \2A 2 k\] 
[(k 2 + £ B ) 2 - 4/l 2 k 2 ] 3 
2 + tj(k) 



A 3 



16kt] 2 (k) VI + t](k) 



(91) 



Defining a new condensate wave function \\f(x) = -\/aA(x), 
we obtain the Gross-Pitaevskii free energy functional 



Tgp[v(x)] 



dx 



vHx) 



72 \ 



dr 2M^ 2M, 



B ' 



jUBlvWI 



Mf{x) 

(92) 
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The coupling g describes the two-body interactions among the 
rashbons. It is given by 



d _ An 2 + tj(k) 
a 1 ~ T[l +T ] (K)] i / 2 ' 



(93) 



For A — * oo, the coupling g goes as g 1.11 /A. Therefore, 
the system is a Bose-Einstein condensate of weakly interact- 
ing rashbons for large values of A/kp. Minimizing the Gross- 
Pitaevskii free energy, we obtain |\|/| 2 = fip/g and the rashbon 
density «b = n/2 — \\\i\ 2 . This is consistent with Eq. d87t . 
which gives n = 2aA 2 . Therefore, the rashbon chemical po- 
tential /^b can be expressed as //b = gn/2 or 



/*b 4[2 + / A 

e F 3n[l+r,(K)fl 2 \k F 



(94) 



At large A/kp, this result is in good agreement with the nu- 
merical results in Fig. [2] For A/kp — > oo, fi B /ep goes as 



B. Superfluid Density 

Now we discuss the superfluid densities, and n'|, which 
are the key quantities to determine the Goldstone mode ve- 
locities. For h = 0, the longitudinal superfluid density equals 
the total fermion density n, n\ — n. However, the transverse 
superfluid density does not. It can be expressed as 



n = n — Ha, 



(95) 



where the spin-orbit coupling induced normal fluid density n,\ 
reads 



Therefore, at large A/kp, the transverse superfluid density is 
suppressed by a factor 2m /M^ < 1. For A — > oo, we have 
Mg/(2m) — > 1.20. This means the transverse superfluid den- 
sity approaches the limit 



nf A 

— -> 0.834 for > oo. 

n kp 



(99) 



In Fig. [3] we show the result of nf/n for various values of 
l/(kpa s ). At unitary, it approaches this limit very fast. 




FIG. 3: The transverse superfluid density nf (divided by n) as a func- 
tion of A/k F for various values of l/(k F a s ). The dashed lines are the 
analytical results 2m/ Mg. 



The physical picture becomes clear if we take a look at the 
phase stiffness J ± and /y. At large spin-orbit coupling, we 
obtain 



Therefore, the transverse superfluid density is always 
smaller than n for A + 0, as shown in Fig. [3] Our result is con- 
sistent with the result of the superfluid density for 3D Rashba 
spin-orbit coupled Fermi gases first reported by Zhou and 
Zhang l35ll . This is in contrast to ordinary s- wave fermionic 
superfluids, where the superfluid density always equals the to- 
tal density at zero temperature. 

To understand why < n, we explore the behavior of nf 
at large spin-orbit coupling. To the leading order of A 2 /|/i| 2 , 
nx can be approximated as 



j - Hi. ~ Hi. j - s - " B 

± Am ~ Mg ' ^ ~ Am~ m\' 



(100) 



This result shows that we recover the phase stiffness for a rash- 
bon superfluid with density «b = n/2 and anisotropic effective 
masses Mg > 2m and M fi = 2m at large A/kp. 

On the other hand, the condensation density hq, which is 
another important quantity for fermionic superfluids, can be 
expressed as 



M 2 
8;r 2 



dk- 



i 



,2A 2 2 



8i 4 ki 



1 

2^ 



(k 2 +£ B )[(k 2 + £B) 2 -4^ 2 ki] 2 ' 



Together with Eq. (81) for n, we obtain 



2m 



2m 



(97) 



(98) 



1 



8 



a\2 ' 



(101) 



At large spin-orbit coupling, we have «o («/2)[l + 
(9(A 2 /|/i| 2 )]. This implies that the condensation fraction 2«o/« 
approaches unity at large A/kp, which is consistent with the 
Bose-Einstein condensation of weakly interacting rashbons. 
In contrast to the superfluid density, we find numerically that 
«o is always an increasing function of A/kp. 
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C. Collective Modes 

For h = 0, the expression for the inverse collective mode 
propagator M(Q) becomes very simple. We have 



1 + 



ra tfi 



k+p^k p 



pa pP 
^k+p^k-pV 



T, 



a/3 



I ca ifi 
tk+n ^k 



'k+p 
pa 
c k+p 

A 2 



k P 



k-p / 

q-a/3 
_/ kq' 



n 



kq' 



a/3 
kq ' 



(102) 



where has a nice property = 5 trj g for /; = 0. The 
expansion parameters A, B, D, R are simplified as 



k+p k-p 



A - — 



1 



Z) = - 



^k 

(g> 2 



y 



/? = — 



vy 

ar=± k 



(103) 



Analytical results can be achieved at large spin-orbit cou- 
pling with the help of these simplified formulas. To the lead- 
ing order of A 2 /\fi\ 2 , the parameters A, B, D, R can be well ap- 
proximated as 



4A 2 d, 



B a 2Aa, D =i d, R^ A-d. 



(104) 



Using the expressions for a and d, we find that at large A, B 2 /A 
goes as B 2 /A ~ A white R goes as R ~ A 2 . Therefore, we have 
B 2 /A » R at large /I and the amplitude-phase mixing term 
dominates the low-energy behavior of the collective modes. 
Using the result g = d/a 2 , we can express the Goldstone mode 
velocities as 



/ g»B 

Ik 

jgns 



ML 




(105) 



These results are just the Goldstone mode velocities of a 
weakly interacting Bose condensate, despite that the bosons 
possess anisotropic effective masses. For A — > oo, the ratio 
cf/c] approaches a universal limit 



-j -> 0.913. 



(106) 



In Fig. |4] we show the numerical results for cf and c|j for 
various values of l/(kpa s ). For A/kp — > oo, they both go as 
1 / yA, independent of the interaction parameter l/(kpa s ). 




FIG. 4: The transverse and longitudinal velocities of the Goldstone 
mode, and c| (divided by the Fermi velocity up = kp/m) as func- 
tions of A/kp. The dashed lines corresponds to the analytical results 
( 1 1 05 b for large spin-orbit coupling. 



Actually, to the leading order of A 2 /|//| 2 , the inverse boson 
propagator M(Q) (at h = 0, M X2 (Q) = M 2 i(0) can be ap- 
proximated as 

Mn(fi) = M 22 (-0 - T\Q) + 2d A 2 , 

M 12 (0 = M 21 (0-dA 2 . (107) 

At large spin-orbit coupling, the two-body vertex function can 
be well approximated as F~'(0 -a[/v„+jUB-WB(q)], where 
WB(q) = qi/(2Mg) + q 2 /(2Mg) is the rashbon dispersion. 
Using the result from the Gross-Pitaevskii free energy, //b = 
gits, the dispersion of the Goldstone boson can be expressed 
as 

co(q) = ^ B (q) [co B (q) + 2gn B ]. (108) 

This is nothing but the Bogoliubov spectrum in a weakly in- 
teracting Bose condensate with anisotropic effective masses. 

On the other hand, the mass gap of the Higgs mode can be 
estimated as 



2a 4;;(*)[1 + t](k)] 2 

Mvt — = A . 

d 2 + t](k) 



(109) 



For A — » oo, Mh goes as Mh ~ 1.04-A 2 . While this may not 
be the true mass gap, we see that in the presence of the spin- 
orbit coupling, the Higgs mode is also pushed up to the large 
characteristic energy scale (~ A 2 ) of the system. 
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IV. RESULTS FOR ZERO ZEEMAN FIELD: 2D CASE 



In two spatial dimensions, the two-body bound state exists 
for arbitrarily weak attraction. Here we show that the spin- 
orbit coupling effect enhances the binding energy. For zero 
Zeeman field, the two-body binding energy 7i B is determined 
by the equation 



z 



k 2 + £ B 



1 



(k 2 + E B ) 2 - 4i 2 k 2 k 2 + e B 



0. 



(HO) 



The solution can also be written as E^/A 2 = 1 + tj(k), where 
k = ln(/ta2D) and t](k) is determined by 



i i yrr^ 

— — arctan — — - In 



7 



(HI) 



The ration Eg / e B can be expressed as 7s B /e B = e 2y+lK {\ +rf)/4. 
The asymptotic behaviors of rj(k) are: (1) tj(k) — » 4e~ 2lc ~ 2y - 1 
for k — > -oo; (2) t](k) — » 7t 2 /(4k 2 ) for k — » +oo. Therefore, we 
have E B — > e B for weak spin-orbit coupling (k — » -oo) and 
Eb 3* e B for strong spin-orbit coupling. 

The effective mass M B of the bound state is given by 



2m \ y (k 2 + 7± B ) 2 + 4i 2 k 2 
Mb/ 2^ [(k 2 + 7± B ) 2 



■ 4^ 2 k 2 ] 



2 12 



?(k 



8^ 4 k 2 



2 + £ B )[(k 2 + £ B ) 2 -4/l 2 k 2 ] 2 ' 



which gives 



2m = 1- (77 -1)7(7?) 

M B 2(7, + 1) [1+7(77)]' 



where 



I{ri) = 



1 In 
2~jjj\2~ 



arctan 



2^ 



(112) 



(113) 



(114) 



For k — > +oo, 7(77) has the asymptotic behavior 7(77) — > /c. 
For /e — » -00, we have 7(77) — » I/77. The ^-dependence of 
the binding energy and the effective mass is shown in Fig. [5] 
Similar to the 3D case, the effective mass M B approaches Am 
for k — > +00. 

The order parameter A and the chemical potential // are ob- 
tained from the gap and number equations. The system can be 
characterized by two dimensionless parameters, ln^pfibD) and 
A/kp, which represent the attractive strength and spin-orbit 
coupling strength, respectively. For numerical calculations, 
it is convenient to employ the following analytical forms for 
the gap and number equations, 



In 



V/J 2 + A 2 - /a 



Jo 



dk 



^ k -A 2 /2) 2 + A 2 ' 



2er 



= ^ 2 + A 2 



+ A\ dk 

Jo 



1 - 



ft - ^ 2 /2 



V(ft - ^/2) 2 + A 2 



(115) 




1.6 





ln(Xa 



FIG. 5: The binding energy E B (divided by e B ) and the effective mass 
Mb (divided by 2m) as functions of ln(/la2D)- 




FIG. 6: The pairing gap A (divided by ep) and the quantity fi + Eb/2 
(divided by e F ) as functions of A/kp for various values of the interac- 
tion parameter ln^pfliD)- The dashed lines correspond to the analyt- 
ical results at strong spin-orbit coupling. 



For A — 0, they give the simple analytical results A = V2e B eF 
and n — e-p - e B /2 J6t] . At large spin-orbit coupling, we have 
7-/ < and A «: which leads to the following analytical 
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results 



where the 2D version of the coupling g = d/a 2 reads 



2na 



— =.-^, (116) 



where the 2D version of a reads 

(k 2 + £ B ) 2 +4i 2 k 2 



1 1 + m 

[(k 2 + E B ) 2 - 4A 2 k 2 ] 2 A 2 4jtt] 



(117) 



Note that in 2D the density n is related to the Fermi momen- 
tum kp by n = k 2 /{2n) = ep/n. The numerical results for A 
and ji are shown in Fig. [6] The results are in good agree- 
ment with the analytical results for large A/kp. However, for 
2D case we find that the quantity fi + E B /2 goes down very 
slowly, unlike the 3D case where it goes as l/A at large A/kp. 
This is because for large A, the rashbon-rashbon coupling g 
goes as g ~ 1/ ln(Aa2D) in 2D rather than g ~ 1 /A as in 3D 
(see below). 

For 2D case, the superfiuid density is isotropic. For h = 0, 
the superfiuid density can also be written as n s = n—nx, where 
the spin-orbit coupling induced normal fluid density nx reads 



At large spin-orbit coupling n,\ can be approximated as 

8i 4 k 2 



(k 2 + £ B )[(k 2 + E B ) 2 - 4A 2 k 2 ] 2 ' 



(119) 



which leads to the result 



n s 2m 
n M B ' 



.1 



Mb 



(120) 



These are just the superfiuid density and phase stiffness for 
2D rashbon condensate. In Fig. [7] we show the result for the 
superfiuid density. At large spin-orbit coupling and/or attrac- 
tion, the result is in good agreement with the above analytical 
result. 

Meanwhile, for large spin-orbit coupling, the expansion pa- 
rameters A, B, D, R can be approximated as 

A =! 4A 2 d, B =! 2Aa, D * d, A 2 d, (121) 

where the 2D version of d reads 

(k 2 + £ B )[(k 2 + E B ) 2 + 12A 2 k 2 ] 



d-2^ 



[(k 2 + £ B ) 2 - 4i 2 k 2 ] 3 



1 2 + (l +77)- 1 +3/(77) 



,! 4 



Snj] 2 



(122) 



At large A, B 2 /A goes as B 2 /A ~ ln(Aa2r>) while R goes as 
R ~ (ln(Aa2D) / A) 2 . Therefore, the amplitude-phase mixing 
term B 2 /A also dominates in 2D. Finally, the Goldstone mode 
velocity c s can be expressed as 



g»B 

M B " 



(123) 



= 2n 



2 + (l + rjy_ +3/(?7) 

[i +im 2 ' 



(124) 



In Fig. [8] we show the numerical result of c s and compare 
it with the analytical result. We note that the coupling g 
is dimensionless and goes as g 6n/ ln(Aa2o) at large A, 
in contrast to the 3D case. The rashbon chemical potential 
jUb = 2ju + E B reads ju B = gn B = gn/2. In Fig. |6] we the 
analytical result for fi B /2 by dashed lines and compare it with 
the quantity ji + E B /2. They are in good agreement at large 
A/kp. The large-/! behavior of the coupling g explains why the 
quantity fi + E B /2 decreases slower than the 3D case. 




FIG. 7: The 2D superfiuid density n s (divided by n) as a function of 
A/kp for various values of ln(£ F a 2 D)- The dashed lines are analytical 
result for large spin-orbit coupling. 




FIG. 8: The 2D Goldstone mode velocity c s (divided by vp = kp/rri) 
as a function of A/kp for various values of ln(£ F a 2 D)- The dashed 
lines are analytical result for large spin-orbit coupling. 
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V. FINITE ZEEMAN FIELD EFFECTS 

Now we turn to the case of nonzero Zeeman field h. In 
the absence of spin-orbit coupling, it is known that the BCS 
superfluidity is completely destroyed at large Zeeman field 
h. In the weak coupling limit, there exists a first order 
phase transition from the BCS state to the normal state at 
hcc — Ao/ V2 where Ao is the gap at h = 0, which is re- 
ferred to as the Chandrasekhar-Clogston limit. Further theo- 
retical studies showed that the inhomogeneous Fulde-Ferrell- 
Larkin-Ovchinnikov state can survive in a narrow window 
between hcc and /ifflo = 0.754Ao. The Zeeman field 
effects on fermionic superfluidity in the whole BCS-BEC 
crossover regime have been experimentally studied in recent 
years lll7ll . Two-component Fermi gases with population im- 
balance (Ni + N\) are realized to simulate the Zeeman field 
effect. Around the unitary point (1 /(kpa s ) = 0), the phase sep- 
aration between the superfluid and the normal phases has been 
observed in accordance with the first order phase transition. 
Despite the rich phase structure in the BCS-BEC crossover, 
one finds that the fermionic superfluidity is completely de- 
stroyed at large enough Zeeman field in the whole BCS-BEC 
crossover regime. 

However, in the presence of spin-orbit coupling, the situa- 
tion is completely different. The crucial observation here is 
that the order parameter A never vanishes even for large h, 
as we will show below. Therefore, the phase structure of the 
spin-orbit coupled system becomes highly nontrivial at large 
Zeeman field h. In the following, we will first study the 2D 
case, since in 2D the results of the bulk superfluid properties 
and the collective modes are most pronouncedly associated 
with the quantum phase transitions. 



At the critical point h = h c , the lower branch E k has a linear 
dispersion near k = 0, i.e., 



E k = v c \k\ + <9(|k| 2 ), k -» 0, 
where the velocity v c can be determined as 

AA 



^Ta 2 " 



(127) 



(128) 



The existence of such a gapless fermionic spectrum causes 
singularities of the thermodynamic functions at the critical 
point Co = 0. 

To be specific, let us consider the thermodynamic potential 
Q(H, h) = Q(n, h, AQx, h)\ where 



Q(/a,A) = ]T U 



+ e B 



+ 6 



(129) 



Note that to obtain the thermodynamic potential Q(yU, h), the 
superfluid order parameter A(//, h), which is regarded as an 
implicit function of fi and h, should be determined by the gap 
equation 



1 + a — 

<Tk/4£? 



1 



k 2 + 



= 0. 



(130) 



We now demonstrate that the infrared singularities caused by 
the gapless fermionic spectrum show up at the fourth deriva- 
tives of the thermodynamic function Q(/i, h) with respect to 
the thermodynamic variables fi and h. To this end, we con- 
sider the following two susceptibilities 



A. 2D System 

1. Quantum Phase Transition 

First, we show that for 2D system, there exist two different 
superfluid phases at h + distinguished by the quantity 



Co = H 



2 + A 2 



(125) 



The superfluid phases with Co > and Co < are linked by a 
quantum phase transition. 

For h + 0, the upper quasiparticle spectrum E k is always 
gapped, while the lower branch E k can have zeros at some 
critical Zeeman field. To show this, we employ the following 
identity 

(E+) 2 (E k ) 2 = (E\ - h 2 - A 2 k 2 f + 4^ 2 k 2 A 2 . (126) 

Since we are considering the superfluid phases with A + for 
A + 0, the only possible zero for E k is located at k = 0. This 
zero appears only when the quantity Co is precisely zero or the 
Zeeman field h equals the critical value h c = ^fi 2 + A 2 . For 
Co + or h + h c , the fermionic excitations are fully gapped. 



d 2 Q(AU0 
dn 2 ' 



Xhh 



<5 2 n0u, h) 
dh 2 ' 



(131) 



Xm is related to the isothermal compressibility and Xhh is me 
spin susceptibility. To obtain their explicit expressions, we 
need the derivatives dA(/j.,h)/dfj. and dA(fj., h)/dh. With the 
help of the gap equation d£2(fi, h, A)/dA = 0, we obtain 



<9A(yU, h) 

dn 
dA(n,h) 

dh 



1 dn(ft,h,A) 
A <9A ' 
1 86n(jj, h, A) 
A dA ' 



(132) 



Here A = d 2 £l(fi, h, A) I OA 2 is one of the expansion parameter 
obtained in Sec. II, n is the total density, and 5n = n<\ — ii[ is 
the spin polarization. We note that the delta function term in 
A vanishes automatically. The explicit expression of n and 5n 
can be evaluated as 



n(ji, h, A) 
5n(ji, h, A) 



Et 



k a=± £ k 



(133) 
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Finally, the two susceptibilities can be evaluated as 
dn(/u, h, A) 



Xfiti 



Xhh 



1 ldn(p.,h,A) 
A \ dA 



d6n(p,h,A) 1 / dSn(/j, h, A) 



dh 



Using the following results 
dEl 

dA 
dEl 

dfx 

d Ik = Jl 

dh ET 
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(134) 
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we find that the expressions of and ^/,/, contain integrals 
of the following type 



J, 



Jo 



(136) 



where 



£ 2 = 



(3 3 = 1 



ft 

ft 
E- 



ft 



(137) 



At h = h c , the integrals Ty (i, j = 1,2,3) are infrared safe 
since the quantities <3, go as A: 2 for k — » 0. Therefore, the two 
susceptibilities and Xhh af e continuous across the critical 
point h = h c . However, the /-th derivative of the suscepti- 
bilities with respect to fi or h contains terms whose infrared 
behavior goes like 



Jr*e £4-2/ r 
« Mk —-l 



dkk 



2-21 



(138) 



For I = 2, infrared divergence shows up, which means that 
the fourth derivative of the thermodynamic potential h) 
becomes divergent at the critical point. Therefore, the third 
derivative of the thermodynamic potential is discontinuous 
across the critical point, where the susceptibilities are continu- 
ous but not smooth. Based on these observations, we conclude 
that the 2D system undergoes a third order quantum phase 
transition at h = h c or Co = even though the superfiuid order 
parameter does not undergo characteristic change. 

On the other hand, we can consider the so-called topolog- 
ical invariant N associated with the fermion Green function 
Q{ito, k). For the present 2D system without time reversal and 
spin rotation symmetries, it is defined as J42* 



N = 



— r- f d 2 kdu> 
in 2 J 



Tr \g 



—Q— 

dk x dk y 



-Tr \Q 



dg- 1 



& 



dky dk x d 



(139) 



Using the explicit form of the fermion Green function, one can 
show that N = for C > and N = 1 for C < 0. Therefore, 
the superfiuid phase with Co < or h > h c is topologically 
nontrivial. The phase transition at Co = or h = h c is a topo- 
logical quantum phase transition. The superfiuid phase with 
Co < can be called a topological superfiuid (TSF), while the 
topologically trivial phase with Co > is a normal superfiuid 
phase (SF). 



< 0.15 



1.2 1.4 




FIG. 9: The pairing gap A (divided by ep) and the chemical potential 
H (divided by e F ) as functions of the Zeeman field h/ep. In this plot, 
we take ln(fcp<J2D) = 2 and A /kp = 0. 5. The dashed line denote the 
critical Zeeman field h c 
h c = 0.5 le F . 



y/fj, 2 + A 2 . The value h c for this case is 



2. Homogeneous System 

Now we turn to study the homogeneous system where the 
total density n = kp/(2n) is fixed. In the numerical calcu- 
lations, we take the attractive coupling as ln(fcpa2D) = 2. 
Change of this value does not lead to qualitatively different 
results. For small enough spin-orbit coupling strength A/kp, 
there exists discontinuity for A and fi, in accordance with the 
first order phase transition and the Chandrasekhar-Clogston 
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FIG. 10: The spin polarization P = (n^ — ti\)/n as a function of h/ep. 



limit in the absence of spin-orbit coupling. However, the pair- 
ing gap A never vanishes once the spin-orbit coupling is turned 
on. Actually, the thermodynamic potential Q(fi, h, A) always 
has a minimum at A + even for large Zeeman field. In this 
section, we are interested in the systems with a large spin- 
orbit coupling strength A ~ 0{kp) which can be realized in the 
cold atom experiments. In Fig. [9] we show a typical numerical 
result for A and ju at A/kp = 0.5. Increasing the spin-orbit cou- 
pling strength does not change our results qualitatively. The 
pairing gap A drops down smoothly with the Zeeman field. At 
large h, we find numerically that A goes as A ~ l/h 2 and n 
becomes largely negative and goes like fi h. 

For the present setting of the attractive strength \n(kpa2o) = 
2 and the spin-orbit coupling strength A/kp = 0.5, we find that 
the topological quantum phase transition occurs at h = h c ^ 
0.5 lep- The chemical potential is smooth across the phase 
transition. It reaches a maximum at h - h c and then drops 
down in the topological superfluid phase. The spin polariza- 
tion §n = «| - «| becomes nonzero once the Zeeman field 
is turned on. In Fig. [TO] we show the the spin polarization 
P = 5n/n as a function of h. We find that P is a smooth func- 
tion of h, in accordance with the fact that the phase transition 
is of third order. We note that the convexities of the A-/j and P- 
h curves change at the critical field h c . In Fig. QTJ we show the 
two susceptibilities Xmi and^;,/, as functions of h/ep. The sus- 
ceptibilities are non-analytical at the critical point h = h c , as 
we expected. Note that the spin susceptibility Xhh is nonzero 
even at h = 0. This is a important spin-orbit-coupling effect 
on fermionic superfluidity/superconductivity l66ll . 

The bulk single-particle excitation gap E g is defined as 

E g = min{£+, E k ] = min{£^}. (140) 

k k 

For h + 0, it does not equal the pairing gap or the superfluid 
order parameter A. In Fig. Q~2] we show the bulk excitation 
gap E g as a function of h. It vanishes only at the critical point 
h = h c . In the SF phase (h < h c ), it is a decreasing function of 
h. In the TSF phase (h > h c ), it is non-monotonic and shows 
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FIG. 11: The susceptibilities^,,, and,^;,;, as functions of h/ep. 
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FIG. 12: The single-particle excitation gap E g = mm^E^, E^} (di- 
vided by e F ) as a function of h/ep. 

a maximum at /; > h c . In Fig. Qj] we show the dispersion of 
the lower fermionic excitation E^ for various values of h/h c . 
One finds that the nodes where the excitation gap E„ opens 
are different for h < h c and h > h c . For h < h c and h > h c , 
the node is located at low momentum, while for h » h c , the 
node moves to k = V2£ F . This can be understood as follows. 
For h » h c , the pairing gap A is very small and the chemical 
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FIG. 13: The dispersion of the lower fermionic excitation E k for 
various values of h/h c . 



It is more interesting to understand the above discussion in 
the helicity representation, i.e., Eqs. (|33]l-(l37]i. For h » h c , 
only the lower band with dispersion £r has a Fermi surface 
k — kp and carries the total fermion density n. In this case, 
the upper band plays essentially no role and only the intra- 
band pairing in the lower band is available. For h — > oo, the 
quasiparticle dispersion is well given by 

E k = yj^k ' K) 2 + |A t (k)| 2 . (147) 

Near the Fermi surface k = kp, it gives exactly Eq. (1144b . This 
means that, at large Zeeman field, the interband pair potential 
A s (k) can be safely dropped and the system can be mapped to 
a spinless p x + ip y superfluid ll67ll where the p-wave pairing at 
the Fermi surface k — kp occurs in the lower band. However, if 
the Zeeman field is not large enough, the mapping to a spinless 
p x + ip y superfluid is no longer legitimate. 



potential becomes negative. To the leading order of A/h and 
A/|/z|, the quasiparticle dispersions can be approximated as 



EJ* A (& + c* % ) 2 +A 2 1 + 



h 2 
x 



(141) 



Therefore, at h » h c , the upper branch has a large gap 
h - fi, and the minimum of the lower branch E k is located 
at the "Fermi surface" k — kp determined by ^ - rj^ — 0. We 
obtain 

~k F = ^2 {a 2 +h + tJa 4 + 2A 2 n (142) 

On the other hand, the total density n can be well approxi- 
mated as 



(143) 



We thus obtain kp V2£ F . Therefore, at h » h c , the gap E g 
of the lower branch is opened at the Fermi surface k = k-p. 
Near the Fermi surface, the dispersion E k can be approxi- 
mated as 



^v 2 {k-h) 2 



+ Ei 



where the Fermi velocity Dp = d£ k /dk\ k= ~ k reads 



vp = V2; 



Vp 



(144) 



(145) 



and the excitation gap E g is given by 



E„ = A 



2A 2 kl 

F 



\2A 2 kl+h 2 ' 



(146) 



3. Collective Modes 

Then we turn to the collective modes. Since the fermionic 
excitations are fully gapped except the critical point h = h c , 
the superfluid density can be expressed as 



Zj /\j 2E? 



IE' 

k k 



24 2 



h 2 E 2 A 2 k 2 ± h 2 A 2 ) E 2 
+ a 1 + -^A + 



^k £ k 



2^ k 



(148) 



Analyzing the infrared behavior of the integrals at h - h c , we 
find that n s is a smooth function across the phase transition. 
In Fig. [14] we show the numerical result for the superfluid 
density. In the SF phase, n s is suppressed by the Zeeman field, 
as we expected. However, in the TSF phase, it turns to be 
enhanced by the Zeeman field! To understand this surprising 
behavior, we first take a look at the limit h » h c . In this 
case, A is very small compared to /; and \p\ and the superfluid 
density can be well approximated as 



= n - — > — 1 + 

Al v 1 d h2 ^ 



©(% - &) 

®(*f - Ik|). 



(149) 



Therefore, for h — * oo, we have n s — > n. As we have shown 
above, for h » h c , the topological superfluid state can be 
mapped to a spinless p x + ip y superfluid where p-wave pair- 
ing occurs at a single Fermi surface k - kp. Therefore, the 
fermion pairing in this phase feels less stress from the Zee- 
man field than in the normal superfluid phase. At large enough 
Zeeman field, since the lower hand carries the total fermion 
density «, we should expect n s — > n. 

Next, the behavior of the expansion parameters A,B,D,R 
can be summarized as follows: (1) A, B and R are continu- 
ous but non-analytical at the critical point h = h c ; (2) D is 
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FIG. 17: The velocity c s of the Goldstone mode (divided by the 
Fermi velocity vp of the non-interacting gas in the absence of spin- 
orbit coupling) as a function of h/ep. 




FIG. 15: The expansion parameters A, B,R, D as functions of /z/ep. 



divergent at /( = h c since one integral in the expression is in- 
frared divergent. The numerical results for these parameters 
are shown in Fig. Q3] The ratio B 2 /(AR) which represents the 



mixing strength between the phase and the amplitude modes 
are shown in Fig. Q~6] We find that the mixing is strongest 
near the critical point, while it becomes negligible for h » h c . 
Since the parameters A, B and R are non-analytical at h = h c , 
the velocity of the Goldstone mode, c s , also behaves non- 
analytically at the critical point. In Fig. [17] we show the 
numerical result of c s as a function of the Zeeman field. It 
is suppressed by the Zeeman field in the SF phase (h < h c ), as 
we expected. However, for h > h c , it goes up rapidly and is 
always enhanced by the Zeeman field. Therefore, the low en- 
ergy collective mode, i.e., the Goldstone sound mode, is a sen- 
sitive probe of the quantum phase transition in 2D spin-orbit 
coupled Fermi gases. For h » h c , the expansion parameters 
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A, B, R can be well approximated as 

a 2 vh i /. if)r 
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(150) 



Since all integrands are peaked at the Fermi surface k — 
kp, we can set k = kp for the integrands except E k 

„2 



^vp(k — kp) 2 + Eg. Meanwhile, we have £k %v.rjv. 



at 



k - kp since A goes as 1/h . Then we obtain 
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Therefore, for ft — > oo, we have B 2 /A -> and R -> 1/(8tt). 
Since the phase-amplitude mixing becomes negligible, the 
Goldstone mode velocity is given by c s = -\/n s /{4R). Com- 
bining the fact n s — > n for /1 — * oo, we obtain 



fp for ft — > oo. 



(152) 



This means that the sound velocity approaches the the Fermi 
velocity vp = kp/m of the non-interacting Fermi gas in the 
absence of spin-orbit coupling and Zeeman field. 

On the other hand, the divergence of D at the quantum 
phase transition indicates that Mh vanishes, that is, the mas- 
sive Higgs mode gets softened near the critical point. In gen- 
eral, the Higgs mode is a resonance. Therefore, we study the 
nature ofits spectral density p(a>,q) = -(l/7r)ImFH(w + ;e,q), 
where the Higgs mode propagator Fe(0 is the massive eigen- 
mode obtained by diagonalizing the matrix M(Q). Here we 
briefly discuss the case of q = 0. From the explicit form 
of the collective mode propagator, we find that the spectral 
density arises when the frequency a> exceeds the threshold 
a>th = 2£g, which is just the two-particle continuum E c (q) = 
mink{£^ +q ^ 2 + E k _ q p} at q = 0. Near the quantum phase tran- 
sition point h = h c , the threshold a)^ tends to zero, while the 
pairing gap A (as well as the critical temperature T c ) remains 
relatively large. This indicates that the Higgs mode spectrum 
arises in the low frequency region near the quantum critical 
point. This may be interesting for future experimental search 
of the Higgs mode in fermionic superfluids. 



B. 3D Case 

The three-dimensional case is different from the 2D case. 
The superfluid state at large h is a gapless superfluid in which 



the lowest quasiparticle has gapless nodes 113011 . This differ- 
ence is due to the existence of the k z degree of freedom. From 
the identity 



(El 



■ ^ 2 ki) 2 + 4A 2 k 2 ± A z , (153) 



we find that the lower branch E k has some zeros located at the 
£ ; -axis (k_L = 0). The k z of these zeros are determined by 



(k 2 /2-vf 



+ A - h 



(154) 



These zeros can be called Fermi points. There exist three 
possible phases according to the number of the Fermi points: 
(1) For h < A, we have a fully gapped phase (SF) without 
Fermi points; (2) For /u > and A < h < tJ/j 2 + A 2 , we 
have a gapless phase (called gSF-I) with four Fermi points 
given by (k ± ,L) = (0,±[2(jj + V/1 2 - A 2 )] l/2 ) and (k ± ,k z ) = 
(0,±[2(ju - V/i 2 - A 2 )] 1/2 ); (3) For h > Vj" 2 + A 2 , we have 
a gapless phase (called gSF-II) with two Fermi points given 
by (k ± ,jy = (0,±[2(/i + V/j 2 - A 2 )] 1/2 ). One can show that 
the quasiparticle dispersion E k is linear around each Fermi 
point, therefore the Fermi point behaves like a Dirac point 
and is topologically protected with a topological charge N a = 
+ 1 i3oTl . The three superfluid phases (SF, gSF-I, gSF-II) are 
therefore separated by two topological quantum phase transi- 
tions. For convenience, in this part we denote the two critical 
Zeeman fields as h c \ — A and h c i — ^/j 2 + A 2 . 

Since the gapless quasiparticle has only Fermi points rather 
than Fermi surfaces, the Dirac delta function terms in the ex- 
pressions of J±,J\\ and A vanish automatically (at zero tem- 
perature). This is quite different from the case of vanishing 
spin-orbit coupling. For that case, the quasiparticle in the gap- 
less phase has gapless Fermi surfaces, and these terms have 
large negative contributions due to the large density of state 
at the gapless Fermi surfaces l52T 460ll . As a result, the gap- 
less superfluid phase in the absence of spin-orbit coupling is 
unstable in a large regime of the BCS-BEC crossover 15311 . 
because A and the superfluid density could be negative. For 
large enough spin-orbit coupling, the gapless phases become 
stable, since these negative contributions are totally removed. 
The cost is that we now have only Fermi points rather than 2D 
Fermi surfaces. 

In Fig. [18] we show the results of the pairing gap A and the 
chemical potential fi for attractive strength l/(kpa s ) = -0.5 
and spin-orbit coupling A/kp = 0.5. Similar to the 2D case, 
we find that the pairing gap never vanishes at large Zeeman 
field, and we do have two critical Zeeman fields h c i 0.34ep 
and h C 2 - 0.47 ep which separate the three superfluid phases. 
For stronger attractive strength and/or spin-orbit coupling, the 
chemical potential fi becomes smaller or even negative. In this 
case, the gapless phase gSF-I with four Fermi points does not 
appear. In Fig. [19] we show the results for the susceptibilities 
and^f/,/,, which are defined in the same way as the 2D case. 
We find that they behave smoothly across the quantum phase 
transitions, which indicates that the quantum phase transition 
in the 3D case is of higher than third order. The reason may be 
understood as follows. Due to the existence of the k z degree 
of freedom, the infrared singularity at kx — * is weakened 
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FIG. 19: The susceptibilities Xmi and^;,;, for the 3D case as functions 
of the Zeeman field. 
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FIG. 21: The transverse superfluid density (divided by n) for the 
3D system as a function of the Zeeman field. 



by the integral over k z . Therefore, the expansion parameters 
A,B,R become smooth at the phase transition points, while 
D does not diverge but exhibit non-analytical behavior. The 
results of these parameters are shown in Fig. [20] At large h, 



we also have B — > and B 2 /A <K R, i.e., the phase mode 
decouples with the amplitude mode at large h. 

Since the gapless quasiparticle has only Fermi points, we 
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FIG. 22: The transverse and longitudinal velocities of the Goldstone 
mode, ct and c!J as functions of the Zeeman field. 



have n s = n for arbitrary h and 



a=± k 

72 



2E2 



a\ 

\2& 



+ h 



2 £4+^2 k 2 A 2 



(155) 



The numerical result for nf is shown in Fig. [21] Its behavior 
is similar to the 2D superfluid density: it decreases at small 
h, and then turns to increase at large Together with the 
behavior of the expansion parameters A, B and R, the results 
of the Goldstone mode velocities cf and c'| are shown in Fig. 
l22l The behavior of the transverse velocity cj is similar to the 
2D case. The longitudinal velocity c| is always an increasing 
function of h, since we have «'| = n for arbitrary h. Unlike the 
2D case, these velocities behaves smoothly across the quan- 
tum phase transitions. 



VI. SUMMARY 

In summary, we have investigated some bulk superfluid 
properties and the properties of the collective modes for at- 
tractive Fermi gases with synthetic Rashba spin-orbit cou- 
pling. Our main results can be summarized as follows: 



(A) For zero Zeeman field, we studied some bulk super- 
fluid properties and the collective modes associated with the 
crossover from the ordinary BCS/BEC superfluid to the Bose- 
Einstein condensation of rashbons. The novel bound state, 
rashbon, exists even for negative s-wave scattering length and 
possesses an effective mass which is generally larger than 
twice of the fermion mass. The behavior of the superfluid den- 
sity and the velocity of the Goldstone sound mode manifest 
the rashbon BEC state at large spin-orbit coupling. Especially, 
we showed that these quantities are universal for A » k F , that 
is, independent of the attractive strength denoted by the pa- 
rameter l/(kpa s ). We also derived the free-energy which de- 
scribes the weakly interacting rashbon condensate and deter- 
mined the rashbon-rashbon coupling g. For large spin-orbit 
coupling A, this coupling g goes as g ~ l/A for the 3D case, 
while it goes as g ~ 1 / ln(/lfl2D) for the 2D case. 

(B) For nonzero Zeeman field, we studied the topological 
quantum phase transitions and properties of the collective 
modes across the phase transitions. For 2D case, we found 
that the second derivatives of the susceptibilities Xw an( l Xhh 
as well as some other thermodynamic quantities are non- 
analytical at the quantum phase transition, which indicates 
that the phase transition is of third order. The singularities of 
the thermodynamic function originate from the infrared diver- 
gence caused by the gapless fermionic spectrum at the quan- 
tum phase transition. As a result, the properties of the col- 
lective modes also show non-analytical behavior at the phase 
transition. The superfluid density as well as the velocity of 
the Goldstone sound mode behave oppositely in the normal 
and topological superfluid phases. In the normal superfluid 
phase, they are suppressed by the Zeeman field as we expect 
from the fact that the Zeeman field serves as a stress on the 
Cooper pairing. However, in the topological superfluid phase, 
they become enhanced by the Zeeman field. Especially, we 
find analytically that n s — > n and c s — > vf in the limit h — > 00. 
This unusual phenomenon may be understood from the fact 
that the system can be mapped to a spinless p x + ip y super- 
fluid of the lower band fermions so that the fermion pairing 
does not feel stress from the Zeeman field. For the 3D case, 
we found that the singularities at the quantum phase transi- 
tions are weakened. However, the behaviors of the transverse 
superfluid density and the transverse velocity of the Goldstone 
sound mode are similar with their 2D counterparts, i.e., they 
are suppressed in the normal superfluid phase while enhanced 
in the gapless superfluid phases by the Zeeman field. 

Spin-orbit coupled atomic Fermi gases have been realized 
at ShanXi University 12(1 ] and at the Massachusetts Institute 
of Technology (MIT) 12111 . by using 40 K atoms and 6 Li atoms, 
respectively. These experiments realized spin-orbit coupled 
Fermi gases with equal Rashba and Dresselhaus spin-orbit 
couplings. While the pure Rashba spin-orbit coupling is hope- 
fully to be realized in the future experiments, it is interesting 
to extend our studies to the present experimental systems. On 
the other hand, the collective mod e sp ectrum can be measured 
by using the Bragg scattering [68| 69] and there are some pro- 
posals for the measurements of the superfluid density lf70ll7lll . 
We hope our theoretical predictions can be tested in future ex- 
periments of spin-orbit coupled Fermi gases. 
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Appendix A: Evaluating the collective mode propagator 

In this Appendix we evaluate the explicit form of M(Q) for arbitrary spin-orbit coupling A and Zeeman field h. To evaluate 
the diagonal element Mn(Q), we decompose the diagonal elements of the fermion propagator as follows 



itit„ — sF" 



iu>„ - sE'l 

s,a=± " k 



icj„ - sE" 



(Al) 



where the quantities C'Jf (k) and C% f (k) are given by 



2sE? 



y— ±a 



sEl - g [{Elf - (f^fYPjjh) - A 2 r y k (-h) 
2sE° (E«f- - (E k a ) 2 



(A2) 



Therefore, Mn(0 and M2i(Q) can be evaluated as 

1 l^^^f(^ a k+p )-f(tEl) 
M„( W , q) = M 22 (-u, q) = T7-^XZ 1 ^ Tr ft** + p)C 22< k " P) • < A3 > 

U Z s,t=±a#=± k W -^k+p +f£ k-p 

Here and in the following p = q/2 for convenience. To calculate the trace in the spin space, we use the following properties of 
the projectors 



Tr[nv^-pW] = Tr K +P (-^<- P (-^)] = \ 

Tr[^ +p (A)yt_ p (-A)] = Tr[^(-A)^ (A)] = 



1 +QjS 



A Z (K - pi) + h 1 



1 +QjS 



^k+p^k-p 
A 2 (ki - P i) - h 2 



^k+p^k-p 



Finally the result can be rearranged as 

Mn(^q) = i + i 2 J] 



aj8=± k 



TV^(k,q) ^(k,q) 
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where the function *W"f (k, q) is defined as 
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The functions ip k and 9?^ are defined as 



1 ±QjS 



^(k* - P i) - / 2 2 



^k+p^k-p 



(A7) 



We can decompose Mn(w,q) as Mn(w,q) = M|j(w, q) + M n («, q), where Mjj(w, q) and M u (cj, q) are even and odd 
functions of at, respectively. Their explicit form read 
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The functions •Wf (k, q) and t/f (k, q) (i =1,2) read 
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For the off-diagonal fermion propagators, we have 
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Therefore Mi2(<2) an d M 2 i(0 can be evaluated as 



W=± t r,/j=± k W ~ i£ k+p + f£ k-p 
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Completing the traces, we find that they can be expressed as 

Mi 2 (w,q) = M^ 2 (a>,q) + /M n (w,q), M 2 i(w,q) = M[ 2 («,q) - iM n (aj,q), 
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Appendix B: Evaluating the Expansion Parameters 

In this Appendix, we evaluate the expansion parameters A, B, D, R and the phase stiffness J ± , J\\ in the zero temperature limit. 
First, A can be evaluated as 
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where the delta function comes from the terms proportional to f(E£ ) - f(E^ ), Taking the derivatives with respect to u, we 
obtain 
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Completing the summation over a,fi — +, we obtain 
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To evaluate the phase stiffness J ± and 7y or the superfluid density and n s , we assume the order parameter takes the 
form ((£>} = Ae 2,qr where A corresponds to the saddle point solution. The superfluid density is defined as the response of the 
thermodynamic potential to an infinitesimal momentum q, that is, 



Q(q) = O(0) + i«tqi + + <?(q 4 )- 



(B4) 



On the other hand, the thermodynamic potential for arbitrary q can be evaluated by making a phase transformation 4> — > e' q r i[r. 
We have 



1 °° 1 
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where 
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Here t,- (i = 1, 2, 3) and to are the Pauli matrices and the identity matrix in the Nambu-Gor'kov space. 

Then we can obtain the explicit form of the phase stiffness using the fermion Green's function Q(ico„, k) at the saddle point 
via the derivative expansion. Obviously, only the I = 1 and 1-2 terms contribute. We note that the I — 1 contribution is just 
identical to the total density n. After some algebra, we obtain 
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The functions X, Y and Z are derived from the Green's function @. X is related to the trace Tr[£?£?]. The result of X reads 
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Here A ± = (iu>„ ± h) 2 - ^ - A 2 - A 2 k]_. Y is related to the trace Tr^cr^cr,] (i = x,y). Due to the angle integration, the nonzero 
contribution reads 

Y(ia> n , k) = {[(/w„ + A + ft)A_ + 2^ 2 ki(^ k + hWm, ~h + ft)A + + 2,* 2 k ± (ft - h)] 
+ [(/w„ +h- ft)A_ - 2i 2 ki(ft - h)][(ia>„ -h- ft)A + - 2A 2 ki(ft + h)] 



+ 2A 2 A + A_ 



[(^„) 2 - ( J B+) 2 ] 2 [(^„) 2 - (£ k ) 2 ] 2 ' 



(B9) 



27 



Z is related to the trace Tr[£?£?cr,£ ; ] (z = x,y). We have 
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The above results can be simplified by using the method of Laurent expansion. We have 
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Completing the Matsubara frequency summation, we obtain 



k o-=± 



k ff=± 
k 2 

2 



»-*ZZ«(|[ + * 



1.2 



2^ 

2 £ k 



2^ 



l-2/(£?) 



2B« 



_L S ech z ^|, 



AT 



IT 



„_V Yk 2 (— sech 2 ^). 



k a=± 

At zero temperature, we have f(E) — > and ^sech 2 ^ — * 6(E). 
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